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' We extensively investigated the critical behavior of mixtures of colloids and polymers via the 

two-component Asakura-Oosawa model and its reduction to a one-component colloidal fluid using 
. accurate theoretical and simulation techniques. In particular the theoretical approach, hierarchical 

reference theory [Adv. Phys. 44, 211 (1995)], incorporates realistically the effects of long-range 
fluctuations on phase separation giving exponents which differ strongly from their mean-field values, 
^v^j . and are in good agreement with those of the three-dimensional Ising model. Computer simulations 

■ combined with finite-size scaling analysis confirm the Ising universality and the accuracy of the 

theory, although some discrepancy in the location of the critical point between one-component and 
full-mixture description remains. To assess the limit of the pair-interaction description, we compare 
one-component and two-component results. 
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Mixtures of colloids and polymers are exciting fluid systems. This is partly due to their industrial importance, but 
also to the fundamental physical insights they provide (Hl^y]- An important mechanism governing the phase behavior 
of colloid-polymer mixtures is the depletion effect, which leads to an effective attraction between the colloids 0, Q 
qq , where the polymer fugacity may be regarded as the analogue of inverse temperature. In particular, for sufficiently 
[~^. ' large size ratio R g /R c , where R g and R c are respectively the radius of gyration of the polymer and the radius of 
I/") , the colloid, a colloid-polymer mixture exhibits stable fluid-fluid phase separation and a solid crystal phase. Since 
colloidal particles can be visualized close to single particle resolution using confocal microscopy, exciting real-space 
investigations of fluid criticality are nowadays possible. In one such experiment, interface fluctuations were directly 
visualized . The snapshots given in Ref. show that the interface fluctuations become more pronounced upon 
approach of the critical point, a consequence of the diverging correlation length. In a more recent experiment, the 
critical exponent of the correlation length was extracted from real-space data and shown to be compatible with the 
three-dimensional (3D) Ising exponent @. This finding is consistent with earlier colloid-polymer experiments @, 0j 
whose results were interpreted in terms of renormalization of the Ising critical exponents [|| • Compared to the wealth 
of results on the overall shape of the phase diagram, little theoretical attention has been given to the critical behavior 
of colloidal dispersions. While the universality class is not expected to be changed from the Ising one, characteristic 
O ■ of simple fluids, the greater flexibility of colloid-colloid interactions with respect to their atomic counterparts may 

allow the study of several nonuniversal aspects of criticality, which are difficult to bring forth in atomic fluids. 
^ . Inspired in part by the above mentioned experiments, the present paper aims to describe fluid criticality in colloid- 
K^J ' polymer mixtures theoretically. In particular we focused on fluid-fluid phase separation and on the critical behavior 
of the Asakura-Oosawa (AO) model 0, 0| via numerical simulation and liquid-state theory. On the simulation side, 
we considered the fluid-fluid phase diagram both of the AO colloid-polymer mixture, and of the one-component fluid 
of particles interacting via the two-body AO potential for a number of size ratios. The critical behavior of the order 
parameter and of the isothermal compressibility were also determined in the binary mixture by finite-size scaling 
techniques. 

Despite its simplicity, extracting the critical behavior of the AO model by theoretical means is challenging. Charac- 
teristic of the mean- field approximation used till now is the parabolic shape of the binodal, corresponding to /3 = 1/2 
where is the critical exponent of the order parameter. Since the AO model belongs to the 3D Ising universality 
class lia, HUH, where [3 « 0.326 0, the true binodal is flatter. In addition, mean-field approximations appreciably 
underestimate the critical polymer fugacity, such that the location of the critical point is not reproduced correctly. 

In this paper we studied the critical behavior of the AO pair potential by the hierarchical reference theory (HRT) 
[l7|. Among liquid-state theories, the HRT has the peculiar feature of implementing the renormalization group 
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machinery, which makes it a reliable tool in the study of the critical behavior of fluids, since it takes into account 
the effect of long range fluctuations on phase separation in a non-perturbative way. This theory is indeed capable 
of getting arbitrarily close to the fluid- fluid critical point. The present analysis complements an earlier investigation 
|18| . where the stress was more on the overall topology of the phase diagram and the stability of fluid-fluid phase 
separation with respect to freezing. Here, we devote our attention solelyto the fluid critical regime. Notice that we 
did not make a direct comparison with the experimental results of refs. , since the mixture considered there has 
a polymer-colloid size ratio larger than one. A one-component treatment where the polymer degrees of freedom are 
traced out is then untenable in that case. 

The main motivations for this paper are the following: first, by comparing the fluid-fluid phase diagram of the 
AO mixture with that of the one-component AO fluid with the effective pair interaction, we aim at assessing the 
accuracy of the latter description as the size ratio is varied. Second, we want to compare the critical behavior of the 
AO mixture as given by finite-size scaling with the results obtained by applying HRT to the AO pair potential. Once 
reduced quantities are used, so as to make allowance for the discrepancy between the critical polymer densities given 
by the one- and the two-component descriptions, we find that the HRT is able to reproduce the asymptotic critical 
power law for the compressibility and the order parameter of the AO mixture with good accuracy. 

The paper is structured as follows: In Sec. IIA the main features of the AO model are recalled. In Sec. HB a 
short overview of the HRT is given. Sec. Ill is devoted to the illustration of the simulation techniques: specifically, 
Sees. IIIA and IIIB deal with the cumulant intersection and the scaling plots methods for determining the critical 
temperature and amplitudes respectively. In Sec. IIIC the unbiased scaling algorithm mentioned above is described. 
A biased version of the algorithm, where one assumes that the universality class in known, is described in Sec. HID. 
In Sec. IV our HRT and simulation results for several size ratios are presented and discussed and our conclusions are 
drawn. 



II. INTRODUCTION TO THE SYSTEM AND THEORETICAL METHOD OF ANALYSIS 



A. The Asakura-Oosawa model 



In this work we considered a binary mixture of colloidal particles and non-adsorbing polymers. Neglecting the 
degrees of freedom of the individual solvent molecules and of the polymer monomers, it is useful to consider effective 
potentials between constituents which are pairwise additive. In this context a very simple and basic model was 
proposed by Asakura and Oosawa 3 (AO model) and independently by Vrij [j| . The colloids are assumed to be hard 
spheres with radius er c /2 and the polymers, with radius of gyration i? g , as inter-penetrating and non-interacting as 
regards to their mutual interactions. However, the polymers are excluded from the colloids to a certain center-of-mass 
distance a cp . Therefore, with respect to their interaction with colloidal particles the polymer molecules are assumed to 
behave as hard spheres of radius er p /2 = i? g , the diameter of the colloid-polymer interaction being a cp = (a c + o~ p )/2. 
The binary AO model is then represented by a binary mixture characterized by the following potentials: 

v pp = 

/ oo r < a cp , , 

oo r < a c 
r > (7 C 

where r is the distance between two particles. The above potentials define what we, in this paper, shall call the 
full mixture description of the AO model, in which the degrees of freedom of both the colloids and the polymers are 
explicitly retained. 

If the polymer degrees of freedom are traced out, the resulting effective (pair) interaction reads as: 



Pvao (r) 



oo r < a, 

3r r 3 



<7 C < r < a c + (Tp (2) 



Q 3 



2(1 + q)a c 2(1 + qfa { 
r > <7 C + a 



here z p is the fugacity of a pure ideal polymer system and q = a p /a c . In Fig. [I]we show the interaction potential 
for several size ratios q. On decreasing q, the attractive well potential becomes deeper and narrower. Here we do not 
intend to describe in detail this model, which is well known in literature (see e.g. |19|). and was previously considered 
also by some of us floj. We just remind that this effective interaction disregards three-body and higher-body terms. 
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FIG. 1: Plots of the effective colloid-colloid pair potential, given by Eq.(J5J, as function of the pair separation r, for several 
colloid-to-polymer size ratios q. 

Interactions beyond the pair level are actually absent for q = ^ < 0.154, but they are present for larger values of q, 
so that vao becomes less and less accurate as q increases. Despite of this limit and of its simplicity, we stress that the 
AO model is able to reproduce also at a fairly quantitative level the main features of the experimental phase diagrams 
for a large range in the size ratios. Finally, we recall that the fluid- fluid dcmixing becomes stable as q increases and 

the depletion potential becomes longer rang ed mmm. 

We applied HRT to the two-body AO potential Eq.(J2J) for different size ratios and we completed the theoretical 
analysis about the critical behavior of this kind of mixtures preliminarily investigated in Ref. [2jj. The theoretical 
results are compared with those we obtained via Monte Carlo simulation and finite size scaling analysis on the full 
mixture as well as on the effective pair interaction Eq.J5J). 

B. Hierarchical Reference Theory 

As introduced before, for the theoretical analysis we consider a model fluid of particles interacting via a pairwise 
additive potential V(ri,r2, . . .r n ) = Yli<j v (\ r i ~ r jl); where is the position of a generic particle i. We assume 
that the two-body potential v(r) is spherically symmetric and results from the sum of a very short-ranged, singular 
repulsion vn(r) and a longer-ranged attraction w(r). In the present case w(r) is «aoM, i.e. Eq.@ for r > <r c . 
The fluid interacting via vpi(r) alone acts as the unperturbed or reference system, whose properties are considered 
as known. Here the reference system is the hard-sphere fluid, whose thermodynamics and correlations are accurately 
described by the Carnahan-Starling equation of state (2^ and the Verlet-Weis parametrization of the two-body radial 
distribution function [24[. We then focus on the perturbation w(r). In HRT 17], w(r) is switched on by taking 
its Fourier transform w(k) and introducing a parameter Q and a potential u>Q(r), such that its Fourier transform 
coincides with w{k) for k > Q, and is identically vanishing for k < Q. As Q evolves from Q = oo to Q = 0, the 
interaction vq(t) = Ufl(r) + wq(t) goes from the reference part vn{r) to the full potential v(r). This procedure for 
turning on the interaction closely resembles that adopted in the momentum-space renormalization group |25l |2(| . At 
a given stage of the evolution, the nature of WQ(r) is such that fluctuations over length scales L larger than 1/Q 
are suppressed, so that critical fluctuations are recovered only in the limit Q — > 0. The corresponding evolution of 
the Helmholtz free energy and n-body direct correlation functions, from the reference fluid to the fully interacting 
one, is described by an exact hierarchy of integral-differential equations. Close to a critical point and at large length 
scales, this hierarchy becomes indeed equivalent to a formulation of the momentum-space renormalization group 27]. 
However, the HRT hierarchy remains valid also away from criticality and over all length scales, thereby describing 
also the nonuniversal behavior of the fluid, which depends on the specific features of the microscopic interaction. 

The first equation of the hierarchy gives the evolution of the Helmholtz free energy Aq of the partially interacting 
system in terms of its two-body direct correlation function in momentum space cq(/c) and the full perturbation w(k): 
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0A Q _ <? J (3) 



dQ 4^ 2 V C Q (Q) 

In the equation above, we have set $(fc) = —w(k)/k-QT, fee being the Boltzmann constant and T the absolute 
temperature, and the quantities Aq, Cq(/c) are linked to Aq and cq(&) by the relations 

Aq = + \ P 2 Mk = 0) 4> Q (fc = 0)] - \ p JJ^L[*(k) - $ Q (k)} (4) 

C Q (fc) = c Q (k) + - $ g (fc) , (5) 

where is the volume and p the number density. These modified free energy and direct correlation function have 
been introduced in order to remove the discontinuities which appear in Aq and cq(/c) at Q = and k = Q respectively 
as a consequence of WQ(k) itself being discontinuous at k = Q. Physically they represent the free energy and direct 
correlation function of the fully interacting fluid as given by a treatment such that the Fourier components of the 
interaction with wavelengths larger than 1/Q are not really disregarded, but instead they are approximately taken 
into account by mean-field theory. In particular, for Q = the modified quantities coincide with the physical ones, 
once the fluctuations have been fully included, while for Q = oo they give the mean-field expressions of the free energy 
and direct correlation function. 

In order to get a closed set of equations, Eq.Q has been supplemented with a closure relation for Cg(fc). This is 
the point where approximations are introduced into the HRT scheme. As in previous applications, our form of CqQc) 
has been inspired by liquid-state theories, and it reads 

C Q (k) = c m (k) + \ Q <5>(k)+G Q (k), (6) 

where chs(^) is the direct correlation function of the hard-sphere fluid, and Aq, GQ{k) are a priori unknown functions 
of the thermodynamic state and of Q. Specifically, the function Qq{k) is determined by the core condition, i.e., the 
requirement that the radial distribution function gc}{r) be vanishing for every Q whenever the interparticle separation 
is less than the hard-sphere diameter a. Aq is adjusted so that Cq(/c) satisfies the compressibility rule. This constraint 
gives the reduced compressibility of the fluid as the structure factor evaluated at zero wavevector, and can be expressed 
in terms of the modified quantities Aq, Cq(&) as 

Cq(* = 0) = ^. (7) 

The compressibility rule Q plays a fundamental role in the present scheme. In fact, when Aq in Eq.@ is determined 
via Eq.lJJJ) and the resulting expression for Cq(/c) is used in Eq.(PJ), one obtains a partial differential equation for Aq 
which reads 



dA Q _cy_ 

dQ 4tt 2 n 



(8) 



where we have set Aq = d 2 AQ/dp 2 , ip(k) = $(k)/§(0) and 

Hk) = CHs(fc) + g Q (k) - [chs(o) + g Q (o)]<p{k) . (9) 

Eq.lJHj is integrated numerically from Q = oo down to Q = 0. At each integration step, GQ^k) is determined by 
the core condition gq(r) = 0, < r < a. This condition acts as an auxiliary equation for determining ip(k) via 
Eq.©. The function Gq(k) has been approximated by a fourth-degree polynomial in the interval < r < a, and 
the equations for the coefficients were obtained. In order to keep the computational scheme relatively simple, further 
approximations were introduced in these auxiliary equations, which amount to decoupling the short- and long-range 
part of the correlations. The details of this procedure have been given elsewhere [28j . 

Since in this work we will be most ly c oncerned with the critical region, it is worthwhile recalling the critical behavior 
of HRT with the closure relation © Hl| ■ In the critical region and at small Q , Eq. (JHJ can be considerably simplified 
in such a way that it depends just on the long-wavelength limit of the direct correlation function. This in turn is 
determined by the compressibility rule Q and the assumption, implicit in Eq.©, that CqIJz) is always an analytic 
function of k: 

C Q (k)~^-bk\ (10) 
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where b is a regular function of Q, p, T. The resulting evolution equation is then cast into universal form by suitably 
rescaling the density and the free energy. The critical behavior of the theory is analyzed in terms of fixed-point 
functions and linearized flow of the rescaled free energy in the neighborhood of the fixed points, along the lines of 
the renormalization group The resulting critical exponents are correct to first order in the expansion in the 

parameter e = 4 — d, d being the dimensionality of the system. In particular, the analytic dependence of the direct 
correlation function on k, also known as the Ornstcin-Zernike ansatz, implies that the critical exponent 77 is zero in 
our approximation. For d — 3 one finds 18] 7 = 1.378, f3 = 0.345, 6 = 5, a = —0.07, v = 7/2 = 0.689, where the 
usual notation for the critical exponents has been used. These exponents satisfy the algebraic relations implied by 
the scaling of the free energy in the critical region [2(|. Below the critical temperature, the theory correctly predicts a 
diverging compressibility inside the coexistence region. However, the compressibility diverges also on the coexistence 
boundary, unlike in the real fluid. This is a consequence of the Ornstein-Zcrnike ansatz (|10fl . 

Finally, we observe that, for the AO fluid considered here, the perturbation $(fc) does not depend on temperature. 
The role of the inverse temperature is instead played by the packing fraction of the polymer in the reservoir 77^, 
defined in terms of the polymer fugacity z p as rf p = irz p <j p /6 (see Eq.@). For 77^ above a critical value rf p cr , the AO 
model will phase separate into a colloid poor phase (the colloidal vapor, with colloid packing fraction ) and colloid 
rich phase (the colloidal liquid, with colloid packing fraction 77'). In this sense, then, rf p is the analogue of inverse 
temperature in fluid- vapor transitions of simple fluids. Both the inverse temperature and rf p appear in Eq.© just as 
external parameters which govern the strength of the interaction for an atomic fluid and the AO fluid respectively. 
Consequently, this does not imply any substantial change in our treatment (see Ref. 20). 



III. SIMULATIONS 



To test the predictions of HRT, large-scale Monte Carlo (MC) simulations of the AO model have also been performed, 
using colloid-to-polymer size ratios q = 0.8, q = 0.56 and q = 0.4. The simulations were carried out in the grand 
canonical ensemble. In this ensemble, the volume V, the colloid fugacity z c , and the polymer fugacity z p are fixed, while 
the number of particles fluctuates. We use cubic simulation boxes with edge L and periodic boundary conditions in all 
d = 3 dimensions. The output of the simulations consists of the distribution Pl(N c \t] p , z c ), defined as the probability 
of observing a system containing iV c colloids at "temperature" 77p, colloid fugacity z c and box size L. 

The simulations of the full AO model, in which both the colloids and the polymers are explicitly retained, were 
performed using the method of Ref. The essential ingredients of this approach are a cluster move [l3l [l4| , a 
reweighting scheme |29j . and histogram extrapolation |30L 13 ll| . In addition, we performed a number of grand canonical 
simulations using the effective pair potential of Eq.(|5J). These simulations, obviously, do not require the cluster move 
of Ref. m and were performed using standard grand canonical MC [32J, |33j . 

One aim of the simulations is to verify to what degree the critical properties of the AO model predicted by HRT are 
reproduced. Previous simulations have shown that the AO model belongs to the 3D Ising universality class [l3l IbH ITU . 
and that pronounced deviations from mean-field behavior become visible upon approach of the critical point. One 
strong point of HRT lies in its ability to yield non-classical (i.e. non mean-field) critical exponents. The latter is 
expected to resemble more closely simulations, and indeed experiments [fJOlQ- In this work, we are particularly 
concerned with the critical behavior of the order parameter and the compressibility in the one-phase region. More 
precisely, defining t = ifplvp cr — 1 as distance from the critical point, the order parameter A = (77J. — rjl)/2 is expected 
to obey A = At 13 , while the (dimensionless) compressibility in the one-phase region \ = v c ((N^) ~ (N c ) 2 )/V is 
expected to diverge as \ — B(-t)^" 1 , with critical exponents (3 and 7, critical amplitudes A and B, and v c = 
the volume of a single colloid. Note that this definition of \ differs from the one adopted in Ref. ITil bv a factor l/v c . 
The quantity x is related to the usual reduced compressibility Xrcd, i-e., the zero-fc value of the structure factor, by 
X = Vc Xicd ■ Preferred values of the critical exponents for the 3D Ising universality class are listed in Table [I] 

To accurately probe the critical properties, the simulation data is analyzed using a number of finite size scaling 
(FSS) techniques, which we will briefly outline in what follows. Most notably, part of our analysis is based on recently 
proposed unbiased scaling algorithms [Tol llll ] . We emphasize here that all FSS algorithms require as input highly 
accurate MC data, typically for a range of system sizes and temperatures. Since our resources are limited, some 



TABLE I: Preferred values of the critical exponents of the specific heat (a), order parameter (/3), compressibility (7), and 
correlation length (v) for the 3D Ising model [l6| . 



a 


P 


7 


z/ 


0.109 


0.326 


1.239 


0.630 
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compromise is unavoidable. The FSS analysis in this work is therefore limited to the full AO model only; no FSS is 
performed on the data obtained using the effective pair potential of Eq. J2J . 

A. Cumulant Intersections (CI) 

The cumulant intersection approach [34[ is a common FSS method to determine the critical temperature nE cr 
from simulation data. Here, the grand canonical distribution PL(N c \rjp, z c ) is used to measure the cumulant ratio 
(to 2 ) /(|to|) 2 with m = N c — (N c ) as function of 7/ p , for a number of system sizes L. In these simulations, the colloid 
fugacity is tuned so as to obey the "equal-area" criterion [3^, |3(| . The data from the different system sizes is expected 
to show a common intersection point, which yields an estimate of the critical polymer reservoir packing fraction rf p cr . 
For the full AO model with q — 0.8, the results of this procedure can be found in Ref. Additional simulations for 
q = 0.56 performed in this work, using three distinct system sizes L — 10, 12, 14 (in units of the colloid diameter), 
show qualitatively similar behavior. The resulting estimates of rf p cr are listed in Table ITT1 Note that the cumulant 
intersection approach can be applied without prior knowledge of the universality class. In this sense, then, it is an 
unbiased method. 

The cumulant intersection method thus requires as input MC data from at least two different system sizes. To verify 
the consistency of the results, however, a higher number is recommended. For q = 0.4, this requirement exceeded the 
computational resources available to us, and so no cumulant intersection result for this size ratio is reported (small 
size ratios q are problematic to simulate due to the high number of polymers these simulations require). The results 
for q = 0.4 are determined by means of the economical scaling algorithm and scaling plots to be described shortly. 

B. Scaling Plots (SP) 

According to finite size scaling theory, the order parameter A obtained in a finite system of linear dimension L 
close to criticality shows a systematic L dependence that can be written as A = L~P/ V M° (tL 1 J v ) , with v the critical 
exponent of the correlation length, and yVf° a scaling function independent of system size |33ll37| . The latter implies 
that plots of L^^A versus tL x l v should collapse onto a single curve. For large tL 1 ^, but still within the critical 
region, this curve should approach the critical power law of the thermodynamic limit. The latter is most conveniently 
visualized using a double logarithmic scale. The data is then expected to approach a straight line, with slope j3 
and intercept equal to the critical amplitude A. This method can thus be used to extract critical amplitudes from 
simulation data. Note that this approach is biased, in the sense that the critical temperature and exponents must 
be known a priori. A possible strategy is to obtain the critical temperature using the cumulant intersection method, 
and to simply assume 3D Ising universality. For many fluids, especially those with short-ranged interactions, the 
latter assumption will be a safe one. Precisely this strategy was followed in Ref. 15 to obtain the critical amplitude 
of the order parameter, as well as the compressibility amplitude, for q = 0.8. The scaling plots obtained by applying 
the same strategy to the q = 0.56 data of this work are qualitatively similar; the resulting estimates of the critical 
amplitudes are listed in Table [D] Note that in these analysis, the colloid fugacity was again chosen to fulfill the 
"equal-area" criterion (see section UlI All . and also that, for the same reason as before, no estimates for q = 0.4 are 
reported. 

C. Unbiased Scaling (US) 

Recently, FSS algorithms were presented that do not require prior knowledge of the critical exponents, nor of the 
critical temperature [ToL ITTl IT3 . l38l l3^| . Instead, these quantities are outputs, and this may prove valuable if there is 
doubt regarding the universality class of a system. In fact, these algorithms were inspired by serious doubts raised over 
the universality class of the restricted primitive electrolyte (which was shown to be that of the 3D Ising model 01). I n 
case of the AO model, there is no doubt regarding the universality class. The motivation for nevertheless using these 
new unbiased FSS methods is a different one. As mentioned in Ref. HoL one problem in simulating asymmetric fluids 
lies in choosing the coexistence chemical potential (or fugacity). A common approach, also adopted by us before, is 
to use the "equal-area" criterion: the fugacity (of the colloids) is chosen such that Pr J (N c \r] I p , z c ) becomes bimodal, 
with two peaks of equal area. Away from the critical point, the peaks in Pi(A^ c |?7 p , z c ) are well separated: equal area 
then corresponds to equal pressure in the two phases 39], a necessary condition for phase coexistence. Close to the 
critical point, however, the peaks in Pl(N c \t] p , z c ) strongly overlap, and a separation in terms of equal area becomes 
rather arbitrary. An alternative procedure is thus desirable. 
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The unbiased algorithm |ToL ITlT ] requires as input the grand canonical distribution Pi,(N c \rf ,z c ) for at least three 
different system sizes L, and for if^ ranging from the non-critical regime toward the critical point. It is therefore 
considerably more expensive than previously discussed methods, which only required data near the critical point. For 
this reason, we consider q = 0.56 only. Starting with 77^ far away from the critical point in the two-phase region, 

the cumulant ratio (to 2 ) 2 / '(to 4 ) is plotted as function of the average colloid packing fraction v c {N c ) /V, with symbols 
defined as before (note that this plot is parameterized by z c ). The resulting curve will reveal two minima, located 
at 77~ and r)J, with respective values Q~ and Q + at the minima. Defining the quantities Q m i n = (Q + + Q~)/2, 
x = Qmin lii(4/eQ m in), and y — (77+ — »7~)/(2A), the points (x,y) from the different system sizes should all collapse 
onto the line y = 1 + x/2. Recall that A is the order parameter in the thermodynamic limit at the considered n^, 
precisely the quantity of interest, which may thus be obtained by fitting until the best collapse onto 1 + x/2 occurs. 
In the next step, r]^ is chosen closer to the critical point, the points (x,y) are calculated as before, but this time A is 
chosen such that the new data set joins smoothly with the previous one, yielding an estimate of the order parameter 
at the new temperature. This procedure is repeated all the way to the critical point, where A vanishes, leading to 
an estimate of 7/ p cr . Moreover, the procedure also yields y as function of x. The latter function is universal within 
a universality class, and for the hard-core square- well (HCSW) fluid can be found in Ref. [n[ Since the AO model 
belongs to the same universality class, we should arrive at a similar result. The latter is verified in Fig. El which shows 
y as function of x obtained in this work, compared to the result of Ref. I 111 Note that, for small x, our data correctly 
approach y = l+x/2. More importantly, the scaling function we obtain agrees well with the one obtained in Ref. I 111 
albeit that our data is "noisier" . This comes as no surprise, since the HCSW fluid is much easier to simulate than 
the AO model, and so higher quality data are more readily generated. The critical value rj^ cr is obtained via the best 
collapse of (t, A) onto a power law. The resulting data is plotted in Fig. H3 (triangles), and the corresponding estimate 
of r/p cr is listed in Tabic [H] The critical exponent (3 and amplitude A are obtained from the slope and intercept of the 
simulation data. For the exponent we find (3 ~ 0.32, which is in good agreement with the accepted 3D Ising value; 
the critical amplitude is listed in Table ITT1 

The unbiased scaling algorithm naturally extends to estimate the coexistence diameter D = (i] l c + i]^)/2 from 
simulation data [Tol Il2|. In contrast to the order parameter, however, the corresponding scaling function for the 
diameter is not universal. In other words, a comparison to Ref. [13 similar in spirit to Fig. is not possible in this 
case. Instead, we simply show in Fig.[3]the output of the scaling algorithm for the q = 0.56 AO model, where we used 
for ?7p cr the value obtained in the previous paragraph. Close to the critical point, the diameter is expected to scale 
as 

D = r) c>CI (1 + aii 2/3 + aat 1 -" + a 3 t) , (11) 

with t the relative distance from the critical point, rj CtCT the critical colloid packing fraction, and non-universal 
amplitudes ai [|lj. The curve in Figgis a fit to the data using the above form yielding ?7 C!Cr = 0.1736, a\ — —0.125, 
<i2 = 1.674 and 03 — —0.875. The HRT diameter (not shown in the figure) varies more slowly with t, its average slope 
being about half of the simulation results in the reduced temperature interval shown in Fig. [21 



D. Economical Scaling (ES) 

The authors of the unbiased scaling algorithm have also presented, in Ref. 0, a biased version of their algorithm. 
If one is prepared to accept 3D Ising universality, the scaling function of the HCSW fluid shown in Fig. |2 which is 
universal within a universality class, can be used to estimate the order parameter of other systems in that class. To 
use this approach, MC data of a single system size obtained close to criticality is in principle sufficient. Since we lack 
the resources to execute the full unbiased scaling algorithm for q — 0.4 and q = 0.8, this economical approach offers 
an attractive alternative. For q — 0.4, a single simulation was thus performed using system size L/a c = 10, while for 
q = 0.8 the data of the largest system of Ref. 13 was used. The triangles in Fig. [7| and Fig. [5] show the corresponding 
order parameter as function of the distance from the critical point. Here, Tf p cr was obtained from the best collapse of 
(i, A) onto a power law; this also yields the critical amplitude A. Using the value of 77^ cr thus obtained for q = 0.4, a 
scaling plot was generated to also extract the critical amplitude B. The latter estimate will not be very precise, but 
seems consistent at least with the trend that smaller size ratios lead to smaller compressibility amplitudes. 

Note that at present, no economical algorithm for the coexistence diameter exists, since the corresponding scaling 
function is not universal in this case Hence, in Table ITT1 no FSS estimate of ?7c.cr could be provided for q = 0.4. 
Instead, we have listed the colloid packing fraction at 7^ cr in the finite system. For q = 0.8, the estimate of Ref. ^] 
is reported, which was derived using the Bruce- Wilding mixed-field scaling method |42fl . 
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0.0 0.1 0.2 0.3 

x 

FIG. 2: Scaling function of the order parameter obtained using the unbiased scaling algorithm. Following the convention of 
Ref. I 111 the scaling function is raised to a negative exponent, with tj> = l/f3. Open circles show results obtained in this work 
for the AO model with q = 0.56; the solid curve is the HCSW result of Ref. El Also shown is the exact small x limiting form 
y= l + x/2. 



IV. RESULTS AND DISCUSSION 

As introduced before, in Fig. ^ we present the pair interaction for the different size ratios we focused on, i.e. 
q = 0.4, 0.56, 0.8. Starting from the effective interaction which describes the AO mixtures, we studied the critical 
behavior by HRT and compared our results with simulations on the same effective interaction. In addition, we 
compared our analysis of the critical behavior with simulations performed on the full, two-component mixture. 

We first considered the whole coexistence region in order to assess the limits entailed by a description of this mixture 
in terms of an effective pair interaction as well as by the particular HRT scheme which we used. In Ref. |2jj, some 
of us showed that HRT yields a very good account of fluid-fluid coexistence curve and of its stability with respect to 
solid phases for several q (q — 0.25, 0.4, 0.6, 0.8). Thanks to its renormalization-group structure, the accuracy of HRT 
in the critical region is remarkable compared to other approximate theories, such as integral equations or perturbative 
methods. In Fig. 0] we show both the data we obtained with HRT and MC simulation on the pair potential vao(?")) 
and the results of our study with MC simulation on the full mixture. There are two main features which emerge 



TABLE II: Summary of our FSS results for the full mixture AO model at various size ratios q. Listed are the critical temperature 
?7p iCr , the critical amplitude A of the order parameter, the critical amplitude B of the compressibility in the one-phase region, 
and the critical colloid packing fraction rj Ct cr . Also indicated is the FSS method that was used to obtain the estimate: cumulant 
intersection (CI), scaling plot (SP), unbiased scaling (US), and economical scaling (ES). The * symbol marks the estimate we 
believe to be the most reliable, in case multiple values are provided. 



<7 




A 


B 




0.8 


0.766 ±0.002 (CI) 
0.7646 ± 0.0003* (ES) 


0.27 ±0.02 (SP) 
0.28 ±0.01* (ES) 


0.053 ± 0.002 (SP) 


0.1340 ± 0.0006° 


0.56 


0.6258 ± 0.0002 (CI) 
0.6256 ±0.0001* (US) 


0.38 ±0.01 (US) 


0.049 ± 0.002 (SP) 


0.1736 ±0.0004 (US) 


0.4 


0.52154 ±0.0001 (ES) 


0.47 ± 0.01 (ES) 


0.045 ± 0.003 (ES±SP) 


~ 0.21 6 



"Taken from Ref. |14| and listed here for completeness. 
'This is no FSS estimate! 
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FIG. 3: Coexistence diameter of the q = 0.56 AO model as function of the relative distance from the critical point, where 
rjp cr = 0.6256 was used. Open circles show simulation results, the closed circle is our estimate of r; CiCr obtained by fitting the 
simulation data to Ea. ljllfl : the curve shows the fit itself. 



from the comparison: on the one hand, the agreement between theory and simulation results on the pair interaction 
is quite good, and increases on increasing the size ratio, because of the theory being tailored, at the present stage, 
to relatively long-ranged interactions (the apparent discrepancy in this trend, central panel, is due to the smaller 
number of MC data close to ?7 ClCr for q — 0.56 than for q = 0.4). On the other hand, the agreement between the 
effective pair interaction description and the full binary mixture increases on decreasing the size ratio, both for the 
overall coexistence curve and the critical packing fractions rf v cr and ?7 c ,cr- This is due to a deficiency of the effective 
pair potential, which is more appropriate to describe colloid-polymer mixtures for small q. In general, the full AO 
mixture appears to be much more sensitive to a change in q than its one-component representation in terms of the 
effective pair interaction «Ao( r )- This is especially true for rjp cr . In fact, most of the deviations of the HRT critical 
parameters from the simulation results of the full mixture are not due to the HRT approximation, but to the modeling 
of the mixture in terms of a pair interaction, as one can appreciate by comparing the HRT coexistence curves with 
the simulation results for the pair interaction model. In Table IHlI the HRT results for the critical packing fractions 
rf v cr and 7y CjCr are compared to the mean- field values. We recall that the mean- field approximation is obtained by 
setting Q = oo in Eq.Q. This corresponds to evaluating the excess free energy with respect to the hard-sphere gas 
a la van der Waals as A — Ahs = \p 2 V Jd 3 rvAo( r )- 

The main purpose of this paper is to study in detail the critical behavior of the AO system with simulations and 
theory and to test the accuracy of HRT in reproducing nontrivial critical exponents. The discrepancy between the 
HRT critical point and the value obtained via finite-size scaling on the AO binary mixture is then at least partially 



TABLE III: Critical temperature 77p, C r 
theory for the different size ratios. 


and critical colloid packing 


fraction ry CiCr determined with HRT and 


mean-field (MF) 




q 


Vp,cr 


*/c,cr 


HRT 


0.8 


0.4825 


0.1895 




0.56 


0.4679 


0.2204 




0.4 


0.4404 


0.2257 


MF 


0.8 


0.4f29 


0.1304 




0.56 


0.4086 


0.1304 




0.4 


0.37f8 


0.1304 
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FIG. 4: Coexistence curve of the AO model for several colloid-to-polymer size ratios q. Closed circles are simulation results 
of the full mixture description, where the crosses mark the location of the critical point obtained using FSS. The triangles are 
simulation results obtained using the effective pair potential of Eq.(2). The open circles are the theoretical predictions obtained 
using HRT. Lines connecting the points serve to guide the eye. The agreement between the pair interaction and the full-mixture 
picture increases on decreasing the size ratio. On the contrary, the accuracy of the HRT with respect to the simulation results, 
both on v(r), increases on increasing q. 



accounted for by adopting the reduced "temperature" variable |?7p/?7p cr — 1|, as customary in the study of the critical 
behavior. Moreover, such a discrepancy is not expected to affect the universal features of the transition. We then 
focus on the critical behavior of the AO model, in particular on the study of the order parameter, compressibility 
and correlation function. In Fig. [3] we present the HRT data for the difference between the reduced densities of the 
colloid on the "liquid" and "vapor" branches of the coexistence curve A = (rj c — rj")/2 as a function of the reduced 
"temperature" rf v /rf p cr — 1 for q = 0.8. A linear fit of our results for log(t) < —4 gives a power-law behavior with 
(3 = 0.37. We recall that the asymptotic value of the exponent (3 predicted by HRT is obtained by linearizing the 
renormalization-group flow induced by the HRT evolution equation (J3J) in the neighborhood of the critical point and 
it is given by (3 — 0.345 01 ■ The theoretical results are compared with the asymptotic Ising-3D power law A ~ At@ 
obtained from finite-size scaling. The critical exponent (3 of the finite-size scaling analysis coincides with the generally 
accepted value (3 = 0.324 (see Sec. EH)- The same figure shows the data obtained by simulating the full mixture. 



11 





FIG. 5: Coexistence curve for q — 0.8 close to the critical point on double logarithmic scales (reduced units). The quantities r/^ 
and r)\. are the colloid packing fraction on the low- and high density branch of the coexistence curve respectively. The triangles 
are data obtained in MC simulations of the full mixture AO model extrapolated to the thermodynamic limit using FSS. The 
dashed line is a fit to the simulation data using the critical power law A = A(rf p /rf pcr — I) 13 , assuming the 3D Ising value 
for (3, and fit parameters ri pCT and A taken from Table UTI Closed circles show the HRT result; the dotted line represents the 
mean-field result. 
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FIG. 6: Coexistence curve for q = 0.56 (reduced units). Symbols defined as in Fig. 



The comparison with the mean-field results (dotted curve) evidences how HRT is able to reproduce the non-trivial 
critical behavior of the mixture. In Fig. Eland Fig.[3we present the same analysis for different size ratios. According 
to both HRT and finite-size scaling, the critical amplitude of the coexistence region increases on decreasing q as a 
consequence of the coexistence curve becoming flatter (see Fig. 0J. However, this trend is stronger for the binary 
system than for the one-component fluid, as can be inferred by comparing the values reported in Table 1TT1 with the 
amplitudes obtained by the linear fit of the HRT results A = 0.41, A = 0.52, A = 0.56 for q = 0.8, q = 0.56, and 
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FIG. 7: Coexistence curve for q = 0.4 close to the critical point (reduced units). Symbols defined as in Fig. 
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FIG. 8: Reduced compressibility Xrcd (upper curves) and correlation length £ (lower curves) for the AO fluid with q = O.i 



q = 0.4 respectively. As a consequence, in the right-hand side of the reduced temperature axis, HRT overestimates 
the finite-size scaling data for q = 0.8, while it underestimates them for q — 0.4. The overall agreement remains 
nevertheless quite good, especially if one considers that the theory contains no free parameters. Far from the critical 
point we recover the mean-field trend. 

In Fig. [SI Fig. and Fig. we consider the behavior of the reduced compressibility Xrod of the colloid on the 
critical isochore in the one-phase region. We report both the HRT results and the asymptotic power law obtained 
from finite-size scaling, x rcd - B'|£|" 7 , 7 = 1.239 O, B' = B/rj c ^ with B given in TableHU We recall the HRT 
result: 7 = 1.378, about 10% larger than the correct one ^j- A linear fit of the data shown in the figures for 
l°g(*) ^ ~4 gives 7 = 1.37. The effect of the interaction range on the critical amplitude of the compressibility is 
weaker than on the order parameter, and is again more evident in the finite-size scaling results for the binary mixture 
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FIG. 9: Reduced compressibility Xrcd (upper curves) and correlation length £ (lower curves) for the AO fluid with q = 0.56. 



than in the one-component fluid. In fact, the data reported in Table ITTl show a weak decrease of the critical amplitude 
on decreasing q, while the amplitude of the HRT compressibility is hardly affected at all in the range of q considered 
here. Therefore, as q is decreased, the finite-size scaling line moves below the HRT points. Also in this case we found 
good overall agreement for each q. In the same figures we have plotted also the HRT results for the correlation length 
£ which governs the decay of the correlations near the critical point. The divergence of £ on the critical isochore 
asymptotically close to the critical point is described by the power law £ ~ C|t| _!/ . The current implementation 
of HRT necessarily gives a vanishing critical exponent rj. Therefore, the critical exponent relation between 7 and v 
becomes just v — j/2. As a general remark, we observe that, while the FSS curves represent the asymptotic power-law 
behavior and therefore give a straight line on the log-log plots of Fig. [SJ-Fig. ^3 the HRT results do not, since in 
the crossover region the corrections to the asymptotic scaling are important. These corrections are nonuniversal, and 
depend on the size ratio q. 

It is worthwhile mentioning that the critical behavior of a colloid-polymer mixture of sterically stabilized silica 
spheres and polydimethylsiloxane dispersed in cyclohexane has been investigated experimentallyhjISt- The isothermal 
compressibility, the correlation length 0, the order parameter, and the interfacial tension |jj of the colloid were 
measured as a function of the polymer concentration, and the critical exponents 7, f3, v were found to be slightly 
higher than the asymptotic Ising values. This was explained in terms of exponent renormalization by the factor 
1/(1 — a), a being the critical exponent of the specific heath [9j. On the other hand, exponent renormalization is not 
found here, where the quantity driving criticality is the reservoir density of the polymer, or equivalently its chemical 
potential. This is a field analogous to temperature in thermal systems, and does not cause renormalization. 

In summary, we have carried out a study of the critical behavior of a simple model of colloid-polymer mixtures by 
both MC simulations and HRT. Critical fluctuations in colloidal systems are expected to have a strong effect on phase 
separation, e.g. lowering the free energy barrier for nucleation of the solid in the fluid phase |44|. More generally, 
there exist many examples which show the importance of a full understanding of the universality class of the systems 
in exam. In a previous paper |15| . one of us studied the interfacial tension of the AO model, which gives indication 
of the strength of capillary waves. The mean-field like behavior and the 3D Ising behavior of the capillarity strength 
are profoundly different; then it is important to reproduce correctly the universality class of the systems in exam. 
The importance of an accurate theoretical approach for this study, where the ability of HRT to describe realistically 
criticality can be very helpful, origins also from the numerical efforts and time cost necessary to the simulation of a 
full colloid-polymer mixture. 
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